Plastid genomes reveal evolutionary shifts in elevational range and flowering time of Osmanthus (Oleaceae)

Abstract Species of Osmanthus are economically important ornamental trees, yet information regarding their plastid genomes (plastomes) have rarely been reported, thus hindering taxonomic and evolutionary studies of this small but enigmatic genus. Here, we performed comparative genomics and evolutionary analyses on plastomes of 16 of the 28 currently accepted species, with 11 plastomes newly sequenced. Phylogenetic studies identified four main lineages within the genus that are here designated the: “Caucasian Osmanthus” (corresponding to O. decorus), “Siphosmanthus” (corresponding to O. sect. Siphosmanthus), “O. serrulatus + O. yunnanensis,” and “Core Osmanthus: (corresponding to O. sect. Osmanthus + O. sect. Linocieroides). Molecular clock analysis suggested that Osmanthus split from its sister clade c. 15.83 Ma. The estimated crown ages of the lineages were the following: genus Osmanthus at 12.66 Ma; “Siphosmanthus” clade at 5.85 Ma; “O. serrulatus + O. yunnanensis” at 4.89 Ma; and “Core Osmanthus: clade at 6.2 Ma. Ancestral state reconstructions and trait mapping showed that ancestors of Osmanthus were spring flowering and originated at lower elevations. Phylogenetic principal component analysis clearly distinguished spring‐flowering species from autumn‐flowering species, suggesting that flowering time differentiation is related to the difference in ecological niches. Nucleotide substitution rates of 80 common genes showed slow evolutionary pace and low nucleotide variations, all genes being subjected to purifying selection.


| INTRODUC TI ON
Attractive fragrance and elegant flowers make species of Osmanthus (Oleaceae) widely used in horticulture, food, and medicine throughout temperate to subtropical regions of the northern hemisphere (Xiang & Liu, 2008). Except for Osmanthus decorus (Boiss. & Balansa), Kasapligil, which is found in the Caucasus Mountains (Green, 1972), the other 27 currently accepted species are found in eastern Asia, from the high-elevation Hengduan Mountains to the low elevations of southeast China, Japan, and Korea (Chang et al., 1996;Green, 1958;POWO, 2021). Osmanthus is characterized by its cymose inflorescences, corolla lobes united in pairs at their base and usually forming a tube, and an androdioecious breeding system (Chang et al., 1996;Duan, Li, Zheng, et al., 2019;Green, 1958Green, , 1972Li et al., 2020).
Genetic studies have discussed the delimitation of Osmanthus using a number of gene fragments (Guo et al., 2011;Li et al., 2020;Lu et al., 2011;Yuan et al., 2010) or a few plastomes Olofsson et al., 2019;Zhao, Yang, et al., 2020;Zhao, Ren, et al., 2020). They found that Osmanthus lies in a generic complex in the subtribe Oleinae (Oleaceae), very close to four genera: Nestegis Raf., Notelaea Vent., Phillyrea L., and Picconia DC. (referred to as NNPP here). Within the genus, only sect. Siphosmanthus and sect. Osmanthus were supported (Guo et al., 2011;Li et al., 2020), although detailed studies on infra-generic relationships, as well as the divergence time of lineages, are still lacking. Therefore, a more widely sampled and robust phylogenetic study is urgently needed to guide downstream evolutionary analysis such as molecular dating and morphological trait evolution. In particular, plastomes have been successfully used to construct the phylogenetic framework of Oleaceae (Dupin et al., 2020;Ha et al., 2018;Olofsson et al., 2019), laying the foundation for our follow-up evolutionary analysis.
Flowering time is considered as a crucial life cycle trait that contributes to fitness in plants (Gaudinier & Blackman, 2020). Changes in flowering time between species may reflect the environmental constraints which limit the direction of evolution in certain species (Wadgymar et al., 2018). The variations in light, temperature, and precipitation caused by changes in elevation are also considered one of the main reasons for the high species richness in eastern Asia (Shimono et al., 2010;Sun et al., 2017). In Osmanthus, significant variations in elevational range and flowering time between species have been noted by previous morphological and ecological surveys . For example, high-elevation species such as O. yunnanensis (Franchet) P.S. Green and O. delavayi Franchet tended to bloom in spring (March to April), while low-elevation species such as O. cooperi Hemsl tended to bloom in autumn (September to October). Exploring the evolutionary patterns of flowering time and elevational shifts within these species can provide insights into how Osmanthus species adapted to different environments (Gaudinier & Blackman, 2020).
Plastomes play an important role in plant phylogeny and evolutionary studies due to them being structurally stable, generally maternal inherited, and with low levels of recombinant DNA (Gu et al., 2018;Wu et al., 2017). Genes of plastomes primarily encode core components of the photosynthetic machinery and factors involved in their expression and assembly, as well as those of self-replication. (Gao et al., 2018;Mohanta et al., 2020;Pilot et al., 2018). Thus, plastomes are generally considered to be conserved, in terms of genomic structures and substitution rates, among the majority of Angiosperms. Recently, however, increasing studies have detected positive selection signals in plastid genes. For example, psaA genes in the 22 closely related Oryza species were found to be undergoing positive selection, which could be related to the adaptation of rice plants to habitats with different light conditions (Gao et al., 2019).
The evolution rates of some genes (psaA, rpl16, ndhA, and ndhH) in Rhodiola plants were also found to be accelerating, which possibly allowed Rhodiola species to adapt to the harsh environment of the Qinghai-Tibet Plateau, such as low carbon dioxide concentration and high solar irradiation (Zhao, Yang, et al., 2020;Zhao, Ren, et al., 2020). Furthermore, the evolutionary rates of matK genes in some low-altitude and recently derived lineages of Dysosma were found to be significantly accelerated, which may reflect the adaptability of these species to new environments (Ye et al., 2018). Thus, genetic content held in plastomes can provide useful information to enhance our understanding regarding adaptive evolution in plants.
Herein, we assembled plastomes of 16 of the 28 Osmanthus species that are currently accepted (POWO, 2021), with 11 plastomes newly sequenced, which allowed us to (1) provide the most well-sampled phylogenetic and molecular clock analyses of Osmanthus to date; (2) study the evolutionary patterns of elevational shifts and flowering time within Osmanthus based on ancestral traits reconstruction; and (3) calculate the substitution rate of plastid genes and explore whether the differentiation of flowering time and elevational range is related to the environmental pressures on plastid genes.

| Plastid genome sequencing, assembly, and annotation
We newly sequenced and assembled the whole plastomes of 11 spe-  (Duan, Li, Zheng, et al., 2019;Olofsson et al., 2019;Zhao, Yang, et al., 2020;Zhao, Ren, et al., 2020). The selection of the outgroups was guided by previous research (Niu et al., 2020;Olofsson et al., 2019), and we The library was sequenced by Nanjing Genepioneer Biotechnologies Inc. (Nanjing, China). Raw reads were obtained and trimmed using CLC Genomics Workbench v9 (CLC Bio, Aarhus, Denmark) with default parameters. The resultant clean reads were then employed to assemble the plastome using the program NOVOPlasty (Dierckxsens et al., 2017) with O. fragrans (GenBank: MG820121) (Duan, Li, Zheng, et al., 2019; as the reference. The resultant genome was annotated by PGA (Qu et al., 2019). The plastomes generated in the present study are available in the NCBI GenBank database. The accession numbers of Osmanthus and outgroups are presented in Table 1.

| Phylogenomic analyses
Two different kinds of plastid data were adopted for the phylogenetic analyses of 16 Osmanthus species and 12 outgroups. The first dataset was made by the concatenation of 80 common protein-coding genes (File S1). The second dataset consisted of whole plastomes, which was designed to be compared with the results of the concatenation sequences. In order to reduce information bias caused by duplicate sequences in the plastomes, one of the two IR regions was removed (File S2). The sequences were aligned using MAFFT (Katoh & Standley, 2013). The poorly aligned sequences were trimmed by trimAL software (Capella-Gutiérrez et al., 2009). The best-fitting nucleotide substitution models for each gene and for the whole plastomes were calculated by ModelFinder (Kalyaanamoorthy et al., 2017) (Files S1 and S3). Phylogenetic trees were constructed by maximum likelihood (ML), Bayesian inference (BI), and maximum parsimony (MP). ML analyses were performed in IQ-TREE v1.6.12 (Nguyen et al., 2015), with 1 million ultrafast bootstrap analyses. BI analyses were conducted using MrBayes v3.2.6 (Ronquist & Huelsenbeck, 2003). The search started from a random tree and used two independent Markov chain Monte Carlo (MCMC) chains run for 1 million generations with sampling of trees every 100 generations. The first 25% of the trees were discarded as burn-in, and the remaining trees were used to generate the 50% majority rule consensus tree. The average standard deviation of split frequencies reached to 0.01, implying convergence of the two runs. MP analyses were performed in MEGA 7 (Kumar et al., 2016) using the search method of tree-bisection and reconnection (TBR) with 1000 random addition replicates. Branch support was calculated based on the bootstrap method with 1000 replications.

| Molecular clock estimation
MCMCtree in the PAML 4.9j package (Yang, 2007) was used to estimate divergence time with an approximate likelihood calculation. The BI phylogenetic tree based on full-length plastomes from 16 Osmanthus and 12 outgroups in subtribe Oleinae was used for divergence time estimation. Estimation of the overall substitution rate was conducted using the BaseML program in PAML. The overall substitution rate (rgene_gamma) and rate drift parameter (sigma2_ gamma) were set as G (1, 128) and G (1, 4.5), respectively. The gradient and the Hessian matrix were estimated using the MCMCtree program in PAML under the GTR substitution model. The MCMC process of MCMCtree was run to sample 20,000 times, with sample frequency set at 500, after a burn-in of 10 million iterations. In total, the MCMC ran for 20,000,000 iterations. Two independent runs were performed to ensure convergence. Distributions of the parameter from MCMC samples were checked by effective sample sizes (ESS) > 200 using Tracer v.1.7 (Rambaut et al., 2018). The results of molecular dating are influenced by many factors, such as the differences between genes, the selection of evolutionary models, and the setting of calibration points. The calibration data for plants are usually derived from limited fossil evidence or secondary calibrations based on previous molecular dating analyses. Therefore, the final results tend to float with the change in calibration schemes (Sauquet et al., 2012). Due to the current lack of fossil evidence in the genus Osmanthus, a secondary calibration approach was adopted. The time calibrations were selected based on previous molecular clock studies of Oleaceae. First, based on age estimations from Besnard et al.

| Elevational range shifts and flowering time evolution
In order to investigate the evolutionary patterns of flowering time and elevational shifts, we reconstructed the ancestral states for 16 Osmanthus and 12 outgroups based on the time tree obtained in the last step. Data on the elevational range and flowering time were obtained based on field observations, literature studies (Chang et al., 1996;Green, 1958Green, , 1972Green, , 2004Li et al., 2020;Xiang & Liu, 2008), and herbarium information stored in the Chinese   We summarized the traits as discrete variables (coded and listed in File S5). For example, the mean elevational range was divided into 500 m and 1500 m; the flowering time was divided into spring (March to April) and autumn (September to October). The estimation of ancestral states was conducted using the R package: Phytools (Revell, 2012). The best-fitting evolutionary models for the discrete characters were selected using the function "fitMk" in Phytools. An equal rates ("ER") model was selected with the highest Akaike information criterion (AIC) weights for both characters. Posterior probability was summarized from 100 stochastic maps and is presented as a pie chart on each internal node of Figure 3.

| Environmental analysis
To investigate the environmental differences between spring-and "extract" function in the R package "raster" (R Core Team, 2021).
The climate information of all sampled points is shown in File S4.
Phylogenetic principal component analysis (pPCA) was employed based on the "BM" method using the "phyl.pca" function in the R package "phytools" to dimensionally reduce the 19 high-resolution environmental variables to two principal component axes, with a consideration of phylogenetic relationships between species, and to detect correlations among them. The data scaling and centralizing were conducted before the pPCA analysis in the R (R Core Team, 2021) platform.

| Substitution rates and nucleotide diversity of plastid genes
To recognize the mutation hotspots for 80 common protein-coding genes across the 16 Osmanthus species, nucleotide diversity (Pi) values were calculated by DnaSP V. 5.10 software (Rozas & Rozas, 1995). To evaluate whether certain genes in specific lineages within Osmanthus (e.g., species with spring-vs. autumn-flowering time or high vs. low elevation) are undergoing positive selection, we performed pairwise Ka/KS calculations on 80 orthologous genes from the 16 Osmanthus species using KaKs_calculator 2.0 (Wang et al., 2010) with the settings genetic code table 11 (bacterial and plant plastid code) and the "YN" method of calculation. If a particular gene in a particular lineage is positively selected for, then we will see that its paired Ka/Ks will be greater than 1, and if it is being subjected to neutral or purifying selection then we will get a value of 1 or less than 1. In our results, some genes had Ka/Ks values equal to "NA," indicating that no synonymous substitutions (Ks) occurred in these genes.

| Phylogenetic analyses
Phylogenetic trees constructed from two datasets (whole plastomes and concatenation of 80 protein-coding genes) using three different tree building methods (BI, ML, and MP) are shown in File S6. The topologies of all trees were highly consistent and similar, but the BI tree constructed from whole plastomes exhibited the highest branch support, and so was used for display and downstream analyses.
In our new phylogeny (Figure 1

| Ancestral state reconstructions
Ancestral state reconstructions are presented in Figure 3a,b.
Flowering in spring was inferred to be the ancestral state in Olea,

| Environmental analysis
In the pPCA for spring-and autumn-flowering species in Osmanthus, the first two axes explained 38% and 31% of total variation, respectively (Figure 3c

| Genome structures and nucleotide variation
The plastome sizes ranged from 155,155 bp to 155,485 bp (Table 1)

| Phylogenetic relationships of Osmanthus
Previous molecular studies have recognized close relationships between Osmanthus and the other four genera in the NNPP clade (Guo et al., 2011;Li et al., 2020;Lu et al., 2011;Yuan et al., 2010). Our plastome phylogeny also confirm their results and suggest a polyphyletic status of traditional Osmanthus. Particularly, our plastome data find the New Caledonian species, O. austrocaledonicus, is in a subclade consisting of Nestegis and Notelaea species. Nuclear-based phylogenetic trees also agree with this finding (Dupin et al., 2020;Li et al., 2020;Olofsson et al., 2019). Morphologically, Asian Osmanthus species are different from species of the NNPP clade and O. austrocaledonicus in their inflorescences being cymose, axillary, with corollas imbricate in bud, as well as wood with a torus (Dute et al., 2010;Li et al., 2020). Therefore, based on molecular, morphological, and biogeographical evidence, we recommend exclusion of O. austrocaledonicus from Osmanthus to ensure the monophyly of the genus.

| Conserved evolution of Osmanthus plastomes
Osmanthus species displayed the same genetic composition and high collinearity, with no recombination or inversion being found. Most of the genetic regions exhibited low nucleotide variability (Pi < 0.004) and almost all the nucleotide variations were non-synonymous substitutions. This indicates that genetic variation within Osmanthus is very limited. When comparing between genes, we found that the ycf1 gene harbored the most nucleotide variability across Osmanthus species. It is worth noting that many recent studies have found ycf1 to be taxonomically informative in phylogenetic research (Dong et al., 2015;Liu et al., 2010), thus making it a good candidate marker as an Osmanthus plant barcode.
By calculating the nucleotide substitution rate, we found the Ka/ Ks values of 80 common genes were less than 1. This indicates that all the genes are undergoing strong purifying selection, which was expected for genes under strong functional constraints like photosynthesis. Purifying selection usually reduces genetic diversity via selective removal of alleles that are deleterious (Cvijović et al., 2018).
More exactly, the functional importance of a protein is a major determinant of its evolutionary rate (Zhang & Yang, 2015). Our Ka/Ks pairwise calculation did not detect the signal of positive gene selection in specific species groups (high elevation vs. low elevation; spring flowering vs. autumn flowering). This indicates that the plastid genes are likely to not be significantly involved in adaptation to altitude or precipitation. Thus, a nuclear genome-wide or transcriptome approach is necessary to study selection pressures on Osmanthus species.

| CON CLUS ION
This is the first well-sampled report on the plastomes of Osmanthus.
Comparative genomic analysis revealed conserved genome structures and low nucleotide polymorphism, although identifying ycf1 as a phylogenetically informative marker that could be used in future studies to gain a more comprehensive overview of Osmanthus ships between species are also much clearer than those resolved in previous studies. Molecular dating of each lineage was also estimated for the first time using two secondary calibration points, with the results suggesting that Osmanthus is a tertiary relict genus.
The origin and diversification of Osmanthus species is speculated to be closely related to the climatic changes and the orogenesis of the Qinghai-Tibet Plateau during the Miocene. Analysis of character evolution revealed the ancestral states of spring-flowering time and low-elevation distribution in the genus. Differentiation of spring-and autumn-flowering time in Osmanthus was also found to be related to the differences in ecological niche between species.
The evolutionary rates of Osmanthus plastomes are very conservative, and so are unable to provide effective genetic information for the study of adaptive evolution. Our results provide a framework for the further study of the systematics and historical biogeography of Osmanthus, including formalizing a subgeneric classification of Osmanthus. A comprehensive molecular sampling of all species should now be focused on in order to achieve these aims.

CO N FLI C T O F I NTE R E S T
None declared. Funding acquisition (supporting); Supervision (supporting).